Connecting microscopic simulations with kinetically constrained models of glasses 
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Kinetically constrained spin models are known to exhibit dynamical behavior mimicking that of 
glass forming systems. They are often understood as coarse-grained models of glass formers, in 
terms of some "mobility" field. The identity of this "mobility" field has remained elusive due to the 
lack of coarse-graining procedures to obtain these models from a more microscopic point of view. 
Here we exhibit a scheme to map the dynamics of a two-dimensional soft disc glass former onto a 
kinetically constrained spin model, providing an attempt at bridging these two approaches. 

PACS numbers: 05.10.-a, 61.20.Gy, 61.43.Fs, 64.70.Pf 



I. INTRODUCTION 

The origin of the onset of ultra-slow dynamics in glassy 
systems, and in particular, glass-forming liquids, remains 
a murky subject, with many competing ideas and tanta- 
lizing clues as to underlying causes, despite years of ef- 
fort by a large community of researchers [l|. Recently 
it has become increasingly clear that dynamical hetero- 
geneities, regions of atypically fast dynamics that are lo- 
calized in space and time, are intimately connected to the 
phenomenon of glassiness 0, & El IS Ifl i becoming increas- 
ingly important at lower temperature scales towards and 
below the glass transition temperature T g 0, H, H, [nj . 

Early ideas about heterogeneous dynamics focused on 
the idea of co-operatively rearranging regions which grow 
with decreasing temperature pd| . Currently, molecular 
dynamics (MD) simulations of supercooled liquids allow 
much greater access to the microscopic details of this het- 
erogeneity [H, HI, [HI • This has included the observation 
of "caging" of particles, and string-like excitations that 
allow particles to escape these cages 0, [lj| [lj| LL21 , which 
has been confirmed in experiments on colloidal glasses 
[ill. However, MD simulations have the drawback that 
it is difficult to reach the low temperatures and long times 
characteristic of the glassy phase. 

An alternative approach to reach low temperatures 
and long times, but without a microscopic foundation, 
is to study simple models of glas siness, kinetically con- 
strained models (KCMs) [H, M EH I2ISS EJ, such 
as the Fredrickson- Andersen (FA) model 20] or the East 
model |25| . or variations such as the North-East model 
[23| , which mimic the constrained dynamics of real glassy 
systems but have trivial thermodynamics. These may be 
viewed as effective models for glasses, in terms of some 
coarse-grained degree of freedom often labelled "spins" , 
also termed a "mobility field" by some authors [23[. In 
Fredrickson and Andersen's original work, they posited 
that the degrees of freedom may be high and low den- 
sity regions, related to earlier suggestions by Angell and 
Rao [201 . Despite these appealing physical pictures, it 
has not been particularly clear to what physical quan- 
tity this "mobility field" corresponds. If KCMs are truly 
effective models of glassy behavior then it should be pos- 
sible to make a connection between some set of degrees of 



freedom, in a MD simulation, for instance, and a KCM. 

In this paper we propose a specific coarse-graining pro- 
cedure to explore whether a link can be made between 
MD simulations and KCMs of glassiness. Previous work 
in this direction found evidence of dynamic facilitation in 
MD simulations [IS H3], however, there was no attempt 
to map the dynamic facilitation onto a KCM. We use 
an approach directly related to the idea of a "mobility" 
field, using the local mean-square displacement (MSD) 
in a suitably defined box to define a spin variable. Re- 
gions with large average MSD correspond to "up" spins 
and those with low average MSD correspond to "down" 
spins. We give specific details of our procedure below. 

We investigate the time and length-scale dependence of 
this coarse-graining. The two characteristic time scales 
are the beta relaxation time scale, tp (corresponding 
physically to the time for relaxation within a cage) and 
the longer alpha relaxation time scale, t a , which corre- 
sponds to the time-scale on which structural relaxation 
of cages occurs. We find that evidence of dynamic facil- 
itation becomes much stronger at longer times of order 
t a , than at earlier times of order tp. We also study the 
effect of changing the size of the coarse-graining box, I, 
in space and consider values 0.02 < l/L < 0.25, where 
L is the system size. With appropriate choices of time 
and lengthscales, we find a clear mapping from our MD 
simulations onto a KCM similar to the 1-spin facilitated 
FA model. This is not what one might naively expect. 
Since we study a fragile glass-former, we expect to find a 
KCM which exhibits super- Arrhenius relaxation, whereas 
the 1-spin facilitated FA model has Arrhenius relaxation 
- the super- Arrhenius growth of timescales is absorbed 
into the coarse graining time, which is most effective at 
capturing kinetically constrained behaviour when it is of 
order t a contrary to expectations based on the idea of a 
mobility field HaHl- 

The demonstration of a coarse-graining procedure to 
translate from a microscopic model to a coarse-grained 
KCM for glasses can shed light on the following: it pro- 
vides a physical interpretation to the "mobility field" ; it 
can give a stronger theoretical justification for the use 
of KCMs to study glassy dynamics; and may open the 
door to further exploration of the link between micro- 
scopic models and long-time features of dynamics, i.e. 
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answering the question: for a given interparticle interac- 
tion potential, how will the dynamics of the glass behave? 

This paper is structured as follows: in Sec. [IT] we gi ye 
details of our molecular dynamics simulations, in Sec. IIIII 
we describe our coarse-graining procedure and present 
the results of coarse-graining MD simulations, and then 
in Sec. |TV]we discuss our results. 



II. MOLECULAR DYNAMICS SIMULATIONS 

We study the dynamics of a well-characterized glass 
former: the binary soft disc model with a potential of 
the form e (^) 12 in two dimensions El El HI- This 
mixture of discs with size ratio 1:1.4 inhibits crystalliza- 
tion upon cooling. We choose a 75:25 ratio of small to 
large particles and cool the system with the density fixed 
at p = 0.85aii. All temperatures and lengths are quoted 
in the standard reduced units of the Lennard- Jones po- 
tential using the small disc diameter cth, as a length 
scale. This mixture has the same glassy characteristics 
as the model previously documented in Ref. (28[. Post- 
equilibration calculations are performed using the re- 
cently introduced iso-configurational (IC) ensemble [29| . 
Our results are qualitatively similar if we follow a single 
trajectory with no IC averaging, but IC averaging gives 
smoother trends as a function of temperature. 




FIG. 1: Averaged mean square displacement from IC simula- 
tions of length T e at T — 0.360. On this time scale there is 
movement throughout the cell. Not all of the the larger groups 
grow into regions of large displacements at longer timescales. 
We show the corresponding spin configuration, with up spins 
represented as black squares. 

In Fig. [TJ we show the MSD for each particle averaged 
over 500 independent trajectories started from the 
same particle configuration in a N — 1600 particle 
system. Each trajectory is evolved for a time r e , where 
r e is the time that it takes for the self-intermediate 
scattering function for the small particles F s (k,t) = 

w (sin (k \n(t) - n(0)\) I (k \n(t) -n(p)\)) Me 



that the form we use has already been averaged over 
angle) to decay by a factor of 1/e - this is roughly 
t a . The k value chosen is that of the first peak of the 
static structure factor. It is clear that there are regions 
of much higher MSD than the average, and that these 
regions are reasonably widely spaced. The important 
question for mapping the dynamics to a KCM is how 
such regions influence the behavior of their neighbors. 

The MSD in the IC ensemble simulations can be seen 
as a measurement of the propensity for a particle to move 
based on the initial configuration. As noted by Widmer- 
Cooper and Harrowell, each trajectory within the ensem- 
ble does not reproduce the same dynamics [29] . The final 
propensity is therefore the composite of a set of trajec- 
tories that is determined solely by the initial particle po- 
sitions. We can follow the change of the propensity in 
time by following a single trajectory and performing IC 
simulations separated by a time r s (which we mostly take 
to be ~ Te). The KCM that we determine is one that is 
obtained from an IC average over 500 trajectories. 



III. COARSE GRAINING PROCEDURE 

There are two parts to our coarse-graining procedure. 
First, we identify the spins that enter in the KCM (using 
the results of the MD simulations). Second we infer the 
dynamics of these spins. Specificly, we construct a model 
which has a Hamiltonian 



i 

where Sj is a "spin" variable on a site i 6 [1, N$] (where 
TVs is the number of sites in the spin model) for which 
up (sj = 1) corresponds to an active region and down 
(sj = —1) corresponds to an inactive region, and J is 
some (yet to be determined) energy scale. The first part 
of the coarse-graining procedure is to find a way to deter- 
mine the separation of regions into up and down spins. 
In general one might also consider terms in the Hamil- 
tonian related to spin-spin interactions (24[, but in their 
simplest forms, KCMs are usually taken to have the form 
in Eq. ([JJ. This implies that at high temperatures there 
are no static correlations, as one expects in a liquid. The 
model Eq. UJ has no interactions and any glassy phe- 
nomenology must come from the dynamical rules that 
govern how spins flip. These dynamical rules are usually 
stated in the form that the probability of a spin flipping 
is dependent on the state of its neighbors [201 - T° be 
more precise, we can note that Glauber rates for flipping 
a spin i are given by El HH : 

ms)=ms){:\\ s f = = _\ (2) 
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which respects detailed balance, and the concentration of 
up spins (with 7if + = 1) is 



and s = (si, . . . , sjv s )- We determine the function fi(s) 
assuming that it has the form f{m), where m is the num- 
ber of up spins on sites neighboring site i, similarly to the 
formulation of the FA model [2(J. We analyze the data 
from the MD simulations to determine /(m). It is desir- 
able that the results be relatively insensitive to the pa- 
rameters entering the coarse-graining procedure, which 
is what we find. 

Our coarse-graining procedure to determine spins and 
sites is as follows: we perform a set of simulations to give 
n ~ 100 timesteps in the IC ensemble. The timesteps are 
chosen to be tg, 0.6r e , and r e to check the coarse-graining 

in time. We find that the fitting form r e = ^^e(^) 
works well over the entire temperature range we consider 
(T = 0.36 to T = 0.96), although for T < 0.48 the form 

/ 1.1 \2 

T e = 0.025e^ T i works equally well. We take each of 
the snapshots of IC averaged particle configurations and 
coarse-grain in space, by dividing the sample into (l/L) 2 
boxes lying on a square lattice [301 ] , so that each particle 
is assigned to a box. We take I to be small and fixed 
to I — 2 most of the time (this appears to be roughly 
the length-scale of the cage-breaking process, and also 
of the order of the dynamic correlation length, £ (28[) 
and assume I to be temperature and time independent. 
Since we have at most 1600 particles in our system, when 
we go to large coarse graining lengths (I > 8), we start 
to get close to the system size and finite size effects are 
important, i.e. l/L > 0.2. 

We associate a spin with each box, either up or down 
depending on whether the MSD per particle in the box 
is larger or smaller than some cutoff. We adjust the cut- 
off so that n j takes its equilibrium value, Eq. ([3]). This 
leaves the freedom to choose the energy scale J. A seem- 
ingly natural energy /temperature scale associated with 
glassy dynamics appears to be that where the relaxation 
time for the small particles starts to stretch more quickly 
and there is a marked onset of dynamic heterogeneity; 
roughly T ~ 0.5. This is also the temperature where the 
scaling between diffusion and relaxation times changes 
[28| . motivating us to choose J = 0.5 (in the same units 
as T). However, we check that our results are robust un- 
der varying this choice (see Fig. [2]). In general we find 
that if J > 0.3 the probabilities of spin flips that we 
identify are identical. 

Now, one of the assumptions that underlies writing 
down Eq. ([T]) is that there are no static correlations. 
Given that the spins we consider are themselves defined 
through dynamics, albeit within a single coarse-graining 
time, it is difficult to be certain that one can extract 
truly static correlations. Nevertheless, wc checked for 
static correlations in the spin model defined above and 



find that at high temperatures there are no static corre- 
lations, whereas for temperatures below about T = 0.45, 
there are "static" correlations on a lengthscale of up to 
two lattice spacings (at T = 0.36), that appears to be 
growing with decreasing temperature. This is in accord 
with the expectation that as the spins diffuse, they lead 
to local relaxations that appear as a static correlation in a 
single snapshot, but can be resolved, for example, by the 
dynamic four-point susceptibility that has been discussed 
extensively in the literature [HI (Hi, [H, H [H [H, H3, H| . 
It is also possible that extra terms involving spin-spin in- 
teractions should perhaps be included in the model as 
these can also lead to static correlations, but the rela- 
tively short correlation length, and the existence of dy- 
namic correlations as discussed above suggests that we 
can ignore these interactions as a first approximation. We 
shall proceed under this assumption of non-interacting 
spins and discuss some consequences of interactions that 
appear to give small corrections to our non-interacting 
results. 
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FIG. 2: P(m) at m — 0, 1, 2, 3, and 4, for spin flips down 
to up as J is varied, with a coarse- graining timescale of r e at 
T = 0.36 with / = 2. 

In Figs. [2] and [3] we show some of our checks on the 
coarse-graining procedure. All of these are at T = 0.36. 
In particular, in Fig. [2] we show how P(m), the prob- 
ability of a spin flip between consecutive time steps 
changes with variations in J at fixed I = 2 and T for 
spin flips from down to up. [At a fixed temperature, 
P(m) oc f{m).) In Fig. [3] we show how P(m) changes 
with variations in I at fixed J and T for spin flips with 
up to down. Comparable results are found for the spin 
flips not shown. Statistical error bars are comparable to 
the size of the symbols. 

We have thus defined our spins. Now we must under- 
stand their dynamics. To do this, we ask the question, 
for an up or down spin with to nearest neighbors that 
are up spins, what is the probability that it will flip in 
a given time-step? This is the way that the classic FA 
model is posed. We display the function /(to) as a func- 
tion of temperature in Figs. 2] and [5] for coarse graining 
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FIG. 3: P(m) as a function of m for spin flips up to down as 
I is varied, with a coarse-graining timescale of r e at T = 0.36 
with J = 0.5. 



times of tp and r e respectively. 

In Figs. Hk) and [5^) we consider m = 0, 1, and 2, 
whilst for clarity, in Figs. HJa) and [5Jd) we show to = 3 
and m — A. The behaviour between these two coarse- 
graining times is quite distinct. For a coarse-graining 
time of tp there is no strong tendency towards kineti- 
cally constrained dynamics, whereas for a coarse graining 
time of r e there are quite strong indications. For a coarse 
graining time of r e , at low temperatures (T < 0.5), /(0), 
/(I ) and /(2) are distinctly different, and there is qual- 
itative agreement between /(to) determined from either 
up to down spin flips or down to up spin flips. Similar 
results are seen for /(3) and /(4). 

Detailed balance implies that /(m) determined from 
either type of spin flip should be the same if the system 
is described precisely by a model of the type in Eq. |T]) . 
In Figs. 2] and [5] there are small but clear differences be- 
tween /(to) determined from the two types of spin flips. 
We believe that there are two sources for this discrepancy. 
Probably most important is the presence of spin-spin in- 
teractions which are not accounted for in Eq. ([T]). Such 
interactions mean that /(to) does not depend solely on 
to, although at low temperatures it is clear that m is 
the most important variable controlling the behaviour 
of /(to). Secondly, the magnitude of the discrepancy 
between rates varies from temperature to temperature 
point. This is likely to be from biases that are forced on 
us by computational restrictions. In principle, one would 
like to equilibrate a large number of independent parti- 
cle configurations, then perform an IC average for each 
initial condition to determine /(to). In practice, it is 
computationally expensive to equilibrate at low temper- 
atures, so we equilibrate one particle configuration and 
then use this to generate subsequent particle configura- 
tions from one member of the IC ensemble. This intro- 
duces a sampling bias that is likely to contribute to the 
discrepancy between the rates. We note that even for a 
single trajectory, there are signatures of kinetically con- 
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strained dynamics, but cleaner results are obtained from 
our IC averaging, and we expect averages over more in- 
dependent particle configurations to yield more precise 
results. Despite these caveats, Fig. [5] clearly illustrates 
that the coarse-graining procedure we have devised gives 
strong evidence of kinetically constrained dynamics, and 
appears to give a mapping of a microscopic simulation to 
a KCM. 

The signatures of kinetically constrained behavior only 
start to show up in the coarse grained spin model for 
T < 0.5. which is around the temperature at which 
the relaxation time for small particles increases quickly 
with decreasing temperature, and is also around where 
the Stokes-Einstein relation breaks down. An important 
point to note about KCMs is that they are expected to 
give the best account of glassy dynamics when nj -C 1. In 
the temperature range we consider, this is only approx- 
imately true, for instance, at T = 0.36, nj ~ 0.2 (with 
J = 0.5, although for larger J, n-f is considerably smaller, 
and the /(to) we have determined are unchanged). The 
requirement of a vanishing number of up spins suggests 
that to extract a specific KCM from our data, one should 
consider the limit that T — ► 0. Figure[5]certainly suggests 
that as T — > 0, /(0) — > 0, but it is harder to determine 
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the fate of /(to) for ra/0. However, our results suggest 
a KCM that might apply in the limit T — > defined by 
(with a constant) 



/(to) 



0, to = 
a, to > 1. 



(4) 



The most important feature of the model that is realised, 
whether it is exactly as in Eq. Q or not, is that it is 
of the one-spin facilitated type [19J. It appears unlikely 
from Fig. [5] that /(l) ^ as T — > as would be required 
for a two-spin facilitated KCM. We have verified with 
Monte Carlo simulations that this model gives timescales 
that diverge in an Arrhenius fashion at low temperatures. 
However, the fact that the timestep in the KCM is also 
strongly temperature dependent (through r e ) leads to a 
fragile behaviour of timescales when expressed in terms 
of the time units of the MD simulations, i.e. r ~ e A l T 
for some constant A. 



IV. DISCUSSION 

There are several goals in trying to map MD simula- 
tions of a glass former onto a KCM. The first is making 
a connection between microscopic particle motions and 
some effective theory of glassy dynamics. A second is to 
determine the long-time dynamics of a given glass for- 
mer at very low temperatures where MD simulations are 
ineffective, say through Monte Carlo simulations of the 
KCM, or in some cases, analytic calculations. 

The mapping of MD simulations to a KCM that we 
have achieved is not what one might naively expect. 
From the discussions in the literature 0, HH , the expec- 
tation would be that for coarse-graining on some small 
lengthscale (as we do), and on a timescale much less 
than t a or r e , one obtains a KCM which has within it 
the physics of the alpha relaxation time. Our attempts 
along these lines are shown in Fig. [H which clearly in- 
dicates that this expected behaviour does not hold. It 
is only when one coarse-grains on timescales of order t a 
that we start to see effective kinetic constraints emerg- 
ing. In order to get fragile glass like behaviour, as seen 
in the MD simulations, from the constant coarse-graining 
time scenario, one would require the KCM obtained from 
the coarse-graining procedure to be multi-spin facilitated, 
since single-spin facilitated models of the type found in 
Eq. (0| are known to have activated dynamics 19]. This 
suggests that there may be alternative coarse graining 
schemes that capture a connection between MD simula- 
tions and KCM. Nevertheless, to the best of our knowl- 
edge, we have demonstrated the first effective mapping of 
MD simulations onto a KCM. The "spins" of our model 
correspond to regions of high or low MSD per particle 
and hence are in the spirit of the "mobility field" [23| ■ 

Given that the mapping we exhibit does not allow us 
to obtain the fragile glass behaviour of the original model 
from our KCM, it is interesting to ask what physics the 
spins that we map to are sensitive to. We believe they 
are sensitive to slow structural relaxation that proceeds 
on time scales even longer than t a . This appears to be 
consistent with recent work by Szamel and Flenner [39j 
that suggests that continuing relaxations beyond t a lead 
to the timescale for the onset of Fickian diffusion to be an 
order of magnitude longer than t a and grow faster than 
t a at low temperatures. The same authors also proposed 
a non-gaussian parameter which differs from that con- 
ventionally used in studies of glass formers This 
non-gaussian parameter has a peak at later times than 
the conventional one, at timescales of order t a , where it 
appears that heterogeneity in the distribution of parti- 
cle displacements is maximal. This might be related to 
why a "spin" definition based on mean square displace- 
ment, as we use here, is most sensitive to physics on 
the timescale of t a . This suggests that in order to make 
a mapping to a KCM more in line with naive expecta- 
tions, one should use quantities that have maximum con- 
trast on timescales which are considerably shorter than 
t a . Such an approach may be a viable way to construct 
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effective theories of glassy dynamics in this and in other 
systems. We have demonstrated a particular numerical 
coarse-graining procedure, and we hope that our results 
may help point the way to analytic approaches to con- 
nect microscopic models of glasses to KCMs, and further 
insight into the glass problem. 
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